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Abstract 

We investigate some aspects of the soliton dynamics in an a-helical protein macromolecule 
within the steric Davydov-Scott model. Our main objective is to elucidate the important 
role of the helical symmetry in the formation, stability and dynamical properties of Davy- 
dov's solitons in an a-helix. We show, analytically and numerically, that the corresponding 
system of nonlinear equations admits several types of stationary soliton solutions and that 
solitons which preserve helical symmetry are dynamically unstable: once formed, they decay 
rapidly when they propagate. On the other hand, the soliton which spontaneously breaks 
the local translational and helical symmetries possess the lowest energy and is a robust lo- 
calized entity. We also demonstrate that this soliton is the result of an hybridization of 
the quasiparticle states from the two lowest degenerate bands and has an inner structure 
which can be described as a modulated multi-hump amplitude distribution of excitations 
on individual spines. The complex and composite structure of the soliton manifests itself 
distinctly when the soliton is moving and some interspine oscillations take place. Such a 
soliton structure and the interspine oscillations have previosly been observed numerically in 
[A.C. Scott, Phys.Rcv. A 26 578 (1982)]. Here we argue that the solitons studied by A. 
Scott, are hybrid solitons and that the oscillations arise due to the helical symmetry of the 
system and result from the motion of the soliton along the a-helix. The frequency of the 
interspine oscillations is shown to be proportional to the soliton velocity. 



1 Introduction 



In the 1970s Davydov 1 proposed a nonlinear mechanism for the storage and transfer of vibra- 
tional energy (intrapeptide vibration Amid-I) in alpha-helical proteins. As a result of the inter- 
action of high-frequency Amid-I vibrations (vibrations of double C-0 bond of peptide groups) 
with the low- frequency acoustic vibrations of the protein, a self-trapping of the Amid-I vibra- 
tion takes place. This idea has attracted a lot of interest, which has increased even further after 
the appearance of a paper [2] in which Davydov and Kislukha demonstrated that the corre- 
sponding system of equations for a molecular chain admits, in the continuum approximation, 
a solitonic solution. This solution describes a self-trapped quasiparticle (a lump of vibrational 
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deformation 0. 

Since then, various properties of such one-dimensional polaron-like self-trapped states have 
been studied in detail both analytically and numerically (see, e.g., [IJ E3 |5] ) . Dynamical prop- 
erties of Davydov solitons and their formation, given various initial conditions of the chain, 
have been investigated in discrete chains and in continuum models. Most of these results have 
been obtained for a single chain. Often they have involved numerical values of the parameters 
that are characteristic of real proteins; thus very often these results have been discussed in the 
context of an a-helix. In reality, however, real a— helices contain three strands, each of which, 
contains periodically placed peptide groups connected by hydrogen bonds. A three-strand model 
for an a-helix was proposed in [7], where the stationary states were first studied. This model 
represents an a-helix as a three-strand structure with three peptide groups per cell in a plane 
perpendicular to the protein axis. Soon afterwards the properties of such soliton states were 
studied analytically in [HI El and numerically in |10( lllj. 

This model does not include a helical structure of proteins, and so afterwards the model was 
improved in |12[ 113] . In ^2] the soliton solutions which do not break the chiral symmetry were 
found analytically. So far, the most complete numerical study of the problem has been presented 
by Scott in |12j . where the formation of a soliton in a linear chain of a finite length had been 
investigated using the initial excitation of a certain form localised on two of the three peptide 
groups at the end of the chain. Scott showed there that, under such conditions, a soliton can 
be formed and that this soliton propagates along the protein with a constant velocity. It has 
turned out that such a soliton has an inner structure and that some interchain oscillations of 
energy take place. These oscillations were compared by Scott with the lines in experimentally 
measured Raman spectra of living cells (see also ^1] ) . Let us add here also that in this numerical 
modelling only one type of initial excitation was used, and, as a result, only one value of the 
velocity of the soliton of a given symmetry was obtained. However, there are several features 
which suggest that this picture is oversimplified. In fact, we expect the dynamics to contain 
some oscillatory features. This is due to the discrete nature of the chain and it can also be 
related to the helical symmetry of the protein and the symmetry of the initial excitation. 

The aim of the present paper is to study the soliton states in an a-helix, to investigate their 
properties and their stability and to look at the dependence of the internal soliton vibrations 
on the velocity of the soliton propagation. In Section 2 a general description of the model is 
given. In Section 3 the elementary excitations of the a— helix are presented. In Section 4 we 
describe the results of our analytical studies of soliton states in the adiabatic approximation 
while the results of our numerical modelling are presented in Section 5. In Section 6 we discuss 
the applicability of the adiabatic approximation and in Conclusions we make further comments 
on the physical relevance of our results. 



2 The general model 

Protein macromolecules are long nonbranched polymer chains which are formed as a result of 
polymerization of aminoacids. Aminoacid residues in such polymer chains are connected by the 
peptide bonds in which four atoms (OCNH) form the peptide group and two a-carbon atoms of the 
residues are placed in one plane. So the backbone of such a polypeptide chain can be described 
as a set of comparably rigid planes divided by methilene groups (-CHR-). Because peptide 
groups (PGs) are bonded with methilene groups by ordinary bonds, a free rotation of PG planes 
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spatial configurations. Thus in particular, it can be rolled into a helix. Such a configuration of 
the polypeptide chain is stabilized by the intrachain hydrogen bonds which are formed between 
a hydrogen atom of a PG and an oxygen atom of the fourth group along the chain. Such a helical 
structure, called a-helix, has 3.6 peptide groups per turn. Thus, the equilibrium positions of 
the repeated units (PG) in an a-helix are determined by the radius-vectors 

ft(o) I ' -> /2vrL _ . ,2ttI.\ al ,,. 

Kr = r e x cos( ) + e„ sin + e z — 1) 

1 V 3.6 ; v v 3.6 ') 3.6 v ; 

where e*j (i = x,y,z) are unit vectors along coordinate axes, a is a period of the helix, r is 
its radius, and I is an integer labeling each group along the polypeptide chain. The nearest 
neighbours (sites / and I ± 1) along the chain are bound by rigid valence bonds and each Z-th 
group in a helix is bound with (I ± 3)-th groups by soft hydrogen bonds forming three spines 
along the helix. 

The three spines along the a-helix are formed by units with numbers: 

l\ =2>n — 1, I2 = 3re, I3 = 2>n + 1, (2) 

or we can write 

I = 3n + (j - 2) (3) 

where j = 1, 2, 3 and n runs from 1 to N with N being the number of PGs in a hydrogen bond 
strand. 

Thus, for the ennumeration of PGs in an a-helix, we can use the two numbers j and n where 
j, a cyclic index modulo 3, indicates the spine of the hydrogen bond, and n ennumerates PG in 
a spine or elementary cells of three PGs from different spines. We can use a different numbering 
of the cyclic index: j' = j — 2 = — 1, 0, +1, or j" = j — 1 = 0, 1, 2, or j = 1, 2, 3. 

Introducing a double index {j, n}, instead of the single number I, the equilibrium positions 
of PGs in an a-helix Q can be rewritten as 

R j,n = r ( ^ cos (— " °i) ~ e v sm (~ " Q i) J + e z(— + Aj), (4) 

where 9j = — 6>o and Aj = ^ — zq. The spines of hydrogen bonds in an a-helix are also 
rolled into a helix of length 5a with 6 PGs per turn. 

Due to the softness of hydrogen bonds, PGs can be displaced and their positions in an a-helix 

are 

where u n are the displacements of the peptide groups from their equilibrium positions (JIJ. 

The potential energy of displacements depends on the distance between the groups and so 
we can perform the approximation of using only the nearest neighbours interaction. The nearest 
neighbours along the polypeptide chain are bound together by rigid valence bonds, much more 
rigid than the hydrogen bond. We can thus assume that the distances between Z-th and (/± l)-th 
groups are fixed while the potential energy of displacements is determinded only by the variation 
of the hydrogen bond length and, in an harmonic approximation, it can be written as 

v = X n/'',,.:,,. ;) - v(Ro)} = £ ^w H (AR jin . jin ^) 2 , (6) 

j,n j,n 



3 



i?o = 4„ a , n -i = l4°n " = V( 2rsin ^/ 6 )) 2 + ( 5fl /6) 2 , (7) 

is the equilibrium length of the hydrogen bond, and 

A z? _ |p £5 I 7? _ (ffi^n ~ Rj,n-l){uj, n ~ %n-l) 

lxJ ^3,n;j,n-X — — -flj,n-l| — -n-0 — 5 (, 8 J 

are its changes due to the small displacements. The total energy of the displacements is the sum 
of the potential energy © and the kinetic energy which is given by the relation 

T = J2 \Mu hn , (9) 

where M is the mass of a PG and i£j jn = —^£- are the velocities of the displacements. 

Due to the assumption that the valence bonds are sufficiently rigid and that the distances 
between the /-th and (I ± l)-th groups are fixed, the three components of the PG displacement 
are not independent. In fact, we have two conditions which correspond to the assumption that 
the distances between each Z-tli PG and its two neighbours, I — 1 and / + 1, are fixed. For small 
displacements this means that the displacement U[ of the Z-th PG is orthogonal to the vectors 
connecting the /-th and the (/ ± l)-th groups: 

ui ■ {Ri - = 0, ui ■ {Ri - R l+ i) = 0. (10) 

_ -Jr) -Jt) 

Let us represent the vector ui using three orthogonal unit vectors & , ej and e z : 

Jr) (r) . Jt) (t) . -> II 

ui = el'u\' +e l 'u\> +e g uf, (11) 
where e z v)j = uf is the longitudinal, (along the a-helix axis) component of the displacement. The 

-» I -(f) (r) Jt) (t) 

transversal component uf = u\ + et'u) is represented through the radial and tangential 
components relative to the axis. Here 



4r) -. ( 2nl . 2vr/ jt) _ . 2nl 2irl 
ej = e x cos + e v sin ), e) = —e x sin + e v cos 



In this case condition (jlOj) takes the form 

— + 2rsin 2 ( — )n, (r) + rsin( — )u, (t) =0, — u, (ll) - 2r sin 2 ( — )u, (r) + r sin( — )u, W = 0. 
3.6 1 v 3.6 ; 1 y 3.6 J 1 ' 3.6 ; ^3.6 ; ' v 3.6 ; 1 

(13) 

From these equations it is easy to find that 

U M = 0, = j , - J. (14) 

1 ' 1 3.6rsin(2^/3.6) 1 V ' 

Thus, there is only one independent degree of freedom of the PG displacement and the vector 
of the displacement be represented as 

^j,n — Sj,nUj t n (15) 

where 



1 



, 2irn 2Ttn \ 2ir 

-a [ sin(- 6j)e x + cos( — 0j)e y J + 3.6r sin(— )e 2 



(16) 
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C=^a2+(3.6sin(^)rj (17) 

is the unit vector which determines the direction of small displacements without changing of the 
valence bond length and Uj tU is the amplitude of the displacements. 

Taking into account this expression, we obtain the following expression for the change of the 
hydrogen bond length (jHJ) 

&Rj,n;j,n-l = l(uj,n ~ Uj,n-l), (18) 

where 

ra ( it 2tt\ 
7= CRA SiU 3 +3S[U 3T6 l (19) 



Thus, the potential energy (JHJ) is 

^ = E \ W i U 3\ri ~ %n-l) 2 , (20) 

where w = ^ 2 wh is an effective elasticity coefficient. 

Therefore, the Hamiltonian of the a-helix vibrations can be rewritten in the form 



Pj> n i 1 I \2 

2M 2 h J ' 



(21) 



where pj tTl are the momentum operators that are canonically conjugate to the operators of the 
PG's displacement Uj >n . 

We now focus on the Hamiltonian for the quasiparticle. The states of the Amid-I vibrations 
of the peptide groups (or extra electron(s)) in the tight binding approximation are described by 
the Hamiltonian 



l \ m / 



(22) 



where / and m run over the 3N values along the polypeptide chain. Here and A[ are, 
respectively, the creation and annihilation operators of the quasiparticle at the Z-th site of the 
chain; L m are the matrix elements of the excitation exchange between sites I and / ± m. The 
matrix elements L m with m being a multiple of 3 describe the energy exchange between the 
PGs of the same spine while the others describe the excitation exchange between the spines. 
For Amid-I excitations in an a-helix the numerical values of L m decrease with increasing m. 
In what follows we will take into account only two the most important terms: L\ = L which 
describes the interspine exchange, and L3 = — J which describes the intraspine one. The signs 
of the corresponding matrix elements are chosen in such a way that they correspond to the 
polypeptide a-helix ^3 |S] 

Using the double index {j,n}: A[ = Aj >n we can rewrite (|22|) as 

H * = E E (^JrAn " JA+ n {A j>n+1 + A hn ^)) + 
n I j 

+ L[AlJA 3 ^ 1 + A 2 ,n)+A+ n (A hn + A 3tn )+A+ n (A 2jn + A hn+1 )]\ (23) 
where n runs from 1 to N and ennumerates the cells on each of the 3 strands. 
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tortion. Due to the softness of the hydrogen bonds and the stiffness of the valence bonds, the 
distance between the n-th and (n ± 3)-th group changes only under distortions of the a-helix. 
So, taking into account the on-site deformation potential only (the dependence of J on the 
distance between the groups is not so essential for an a-helix 0), we can write the interaction 
Hamiltonian in the form 

Hint = X(Uj,n+l ~ U jtn -i)A^ n A jtn . (24) 

where x 1S a constant parametrising the strength of the exciton(electron)-phonon interaction. 
The total Hamiltonian 

H = H e + H v + H int (25) 

where H e , H v , and H; mt are given by ([2*H|) . (|2*Tj) . and (|24[). respectively, describes the dynamics 
of a molecular helical chain in which the equilibrium positions of its units (PGs) are given by 
the radius-vectors 

= r fe x cos(^j) + e y sin(-^j)J + e z a(n + |) (26) 

where j is a cyclic index (of modulus 3), which ennumerates the three spines along the z axes, and 
n is the position index of an elementary cell within the three strands. Equation (j26j) describes 
a molecular chain which is rolled into a helix with three units per turn of the helix. We can 
consider such a molecular chain as a model of the a-helical protein. In this case we can neglect 
the rolling of the hydrogen bonds into a superhelix. 



3 Elementary excitations in an a-helix 



The quasiparticle Hamiltonian (|23|) can be diagonalized by the following unitary transformation 

= ^= E e lkm v h n(k)B,, k , v h n(k) = (27) 

where the wave number k and the band index ^ are given by 

2vr N-l 2tt 

k = — I, Z = 0,±1,...,±— ^— , M=y^, i/ = 0,±l. (28) 



Under transformation (|27|). the Hamiltonian (|23ft transforms into 



H e = J2Mk)B+ k B„ tk (29) 



where the energy dispersion in the three bands (fi = 0, ±3-) is given by 



^(fc) = J B -2Jcos(A;) + 2Lcos(^+ / uj (30) 



or, in explicit form, by 



E (k) = E - 2Jcos(k) + 2Lcos (^j , 

E±(k) = E - 2Jcos(k) - Lcos Q") ± v^Lsin (31) 



_1_ _■_ XI 111111 \ S 111 (111 J. J. ^ U.VUVJ1 1 V_<U 1 1 1 V t V ^.y V'll V.l V- ll V., 1 UUVlllUlUlVllU W 1 U llv 1 V_J 111 U11\J kj J_/ lllV'IJ Wl 1 1 

bonds in an a-helix. 



Next, we perform a unitary transformation of the lattice variables: 



ft» = -^E^(^) 5 K.-4-,). 

re the 

wavenumber g and frequency 



where aj _ ? and aj :Q are the operators of creation and annihilation of acoustic phonons with 



Uq = 2v a \sm-\, "a = J^- (34) 



As it is convenient to describe the lattice oscillations in the helical symmetry representation 
that we have introduced for the description of excitons, we define the operators b vq as 

a n = E v iA ( l) h ^ ( 35 ) 

V 

where the Vj u (q) were given in the description of the excitons (|27j) . In this formulation, the 
displacement operator is given by 

i 

2 



qv 

and the Hamiltonian H v (|21j) takes the form 

H v = Y,^M, q Ki + h- (37) 



vq 



Thus, the elementary excitations are given by the phonons which correspond to the deformational 
oscillations of the lattice, and the excitons which describe the internal Amid excitations of the 
PG. As an elementary cell contains 3 PGs, the spectrum consists of three exciton bands which 
correspond to the Davydov splitting. This band structure is shown in Fig. 1 for L=12.4 cm -1 
and J=7.8cm™ 1 , which correspond to the a-helix values. 

Finally, we rewrite the interaction Hamiltonian Hi n t (|24l) as 

H int = 7= E {x(q)Bl + ^ k+q B^ k b u , g + X*(g)B£ fc iW+A.} (38) 

* kq^iv 

where 

X{q) = lX \Jlk) 1 2 S[n{q) - (39) 
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Figure 1: The three energy bands Q31|) for J = 7.8cm 1 and L = 12.4cm 

4 Equations in the adiabatic approximation 



In the adiabatic approximation the wavefunction of the system with one quasiparticle is repre- 
sented as 

m)) = u(t)\Mt)), (40) 

where U (t) is the unitary operator of the coherent molecule displacements 



U(t) = exp 



v,q 



(41) 



|l&e(t)) = E^fc(tK,fc|0). (42) 
with functions ipu k(t) that satisfy the normalisation condition: 

£hM*)l a = 1 - (43) 

The coefficients /3uq(t) in (|41jl are, at this stage, arbitrary functions which will be determined 
below. 

In the adiabatic approximation the equations for ^fc(t) and (3 u ,q(t) can be obtained either 
directly from the time dependant Schrodinger equation or as Hamilton equations for the gen- 
eralized variables tp^kit), Pu,q(t) and their canonically conjugated momenta (—i/f^ipl^it) and 
{—i/fyPt qif) by considering 

n = (mm = E MWiki>*k + E + 1) 

+ 7^ E x(g)^,*+ 9 ^fcC^,9 + ^ 1 - g ) (44) 
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= E,{k)^ k {t) + Y, 2 -^=Q^t)^ k -M (45) 

ih ^g = ft^p^ + _L= x*(q)^ tk ' t P^+u,k+q- (46) 
ai VoiV m 

In the first equation, Q u (q,t) is given by 

Qufat) = (^A 2 (#,,, + 0Z u ,- g ) ■ (47) 

In fact, the equation for f3 v>q (t) becomes the equation for Q u (q,t) and takes the form: 
d 2 Q u (q,t) 2 2ixsin<2^ 

-p + WgQ w (g, t) = 'P^kim^k+iW- (48) 

In these expressions, and in what follows, the index v labeling tp and Q is defined modulo 3. 

Next, we seek the stationary solutions of these equations by requiring that 

VVfcW = e-^ G W +fa 'W) (49) 

This immediately tells us that 

Q v {k,t) = e~ ikz ^ Q v (k). (50) 
Here the parameter z s (t) corresponds to the centre of mass of the excitatiton. 

Substituting this ansatz into our equations, we obtain 

[Ml + HVk - E,{k)} %{k) = £ ^Sq^^ - g), (51) 

(u, ? 2 - V\ 2 )QM, t) = £ r,m, +v (k + g), (52) 

where = and = ^f? is the velocity of the propagation of the excitation measured in units 
of the lattice constants. 

Taking into account Q52|) . we see that (|51|) is a nonlinear integral equations. From Q51|) we see 
that tp^,(k) has a maximum at the carrying wave number fc c which corresponds to the minimum 
of HO, + HVk — En(k), i.e. k C n and the excitation velocity V are connected by the relation 

Next, we assume that, in the space representation, the solution is given by a wave packet 
broad enough so that it is sufficiently narrow in the k representation. This means that ip^,(k) 
are essentially nonzero only in a small region of values of k in the vicinity of fc c „. In this case 
we can use the following approximation 



m + nvk - EJk) = A- n( ' k n kc ^ , (54) 

2m M 
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h 2 d 2 E (k) 

A = [Ml + hVk - E,{k)] k=kcii , — = -^PU= fecM - (55) 

To solve 1)51(1 . we introduce the position dependent functions 

^) = ^E ei(fc ~ fccM NVW- (56) 

Note that at x = n this is a unitary transformation of ip^(k) to the site representation. Using 
approximation (|54[). one can transform ()51() into a differential equation for y M (x): 

A^(s) + " E V.CxJe-^**-^))*^^) = 0, (57) 

where 

V v (x) = -j=2L= ^^>Q v (q)mnq. (58) 

Note also that (|57j) is only a zero-order approximation of 1)51(1 and so it corresponds to the 
continuum approximation. Only in this approximation the soliton velocity V and frequency ft 
are constant and the soliton centre of mass evolves with time as z s (t) = Vt + zq and 0(t) = fit. 

Finally note that, when transforming (|51|) into 1)57(1 . one has to be careful with the double 
summation (J2k q ■••)• The wavenumbers k and q are in the first Brillouin zone, —it < k ,q < it, 
and the wave number k — q also has to be in this zone. This is the case for small values of k and q 
(normal processes in exciton-phonon interactions). However, when k and q are close to the edge 
of the first Brillouin zone, it is possible that \k — q\ > it (Umklapp processes). In this case it is 
nessesary to reduce the wavenumber k — q to the first Brillouin zone using the reciprocal lattice 
wavenumber g = 2ir. This does not change the discrete equations due to the periodicity of the 
functions in the space of reciprocal lattice vectors, but is essential when introducing continuous 
functions for the analytical investigations. The Umklapp processes lead to the appearence of 
additional terms in (|57|) for which the double summations are performed in the regions near 
the edges of the Brillouin zone where \q — k\ > ir. The assumption that tp^ and Q(q) are small 
in these regions allows us to consider these terms as a perturbation. Here we do not take this 
perturbation into consideration. A detailed analysis can be found in ^S] where it has been 
shown that allowing Umklapp processes in k space leads to the appearence of a periodical (with 
a period of a lattice constant) Peierls-Nabarro potential barrier for the motion of the soliton 
centre of mass ( ^3 El E2 ) • As a result, in discrete lattices, the "instantaneous" soliton 
velocity depends on time and has an oscillatory component with a period 

2?r . . 

T d = (59) 

"av 

where V av is the average velocity of the soliton propagation in the chain. 

Having found the solutions of (|57|). we can then use the transformation (|27|1 and, taking 
into account (|49|) and (|56|). write down the probability amplitudes for the distribution of the 
excitations in an a-helix: 

% n (t) = ^E e *X^)^(*) = ^X) e ^ ( ^ V ^ )t+ ^ n+<(AH "^ ) VA.(n + ^'- Vt-zo). 
v N a* 6 

(60) 
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Q ° {q) = K -VW3MN g ****** ( 61 ) 

Q+(q) = ^ _ 2 y 2 ^ 3MN E (r (k)Mk + q)+ r+W-(k + q)+ ip*_{k)M k + ?)) > 

(62) 

Q-(<z) = <&(-?)• (63) 
Substituting these expressions into ()58[). we obtain the potentials V^(x). For example, 

= E 3M S S -yl 2 ) e "^^ (fc)e ' (fc+9)x ^ (fc + q) - (64) 



We can see from (|61j) that Qo(q) is essentially nonzero only at small values of q. So we can 
use the long-wave approximation and write the phonon dispersion relation (|H4jl as uj q ~ i/ <x |q r | . 
In this case Vq(x) is given by 

Vo(x) = -Mi^)^^ )|2 (65) 

where v = V/v a is a velocity in units of sound velocity v a . Similarly, the potentials V±(x) are 
quadratic in <pa(x). Therefore, the system of equations (|57jl is a system of nonlinear Schrodinger 
equations (NLSEs). 

We observe that equations (|57j) admit three types of ground state solutions of a soliton type 
which preserve the helical symmetry of the system. Such solutions describe solitons which are 
formed by excitons from only one of the three excitonic bands, i.e. only one function 93^ 7^ 
for a given fi is nonzero and the other two ip v = with v 7^ \i. In such states, according to 
(|61|) -([62 |) . only the total symmetrical distortion of the a-helix takes place, i.e. Qo(q) 7^ and 
Q±(q) = 0. Taking into account Q65|>. we note that these types of solitons are described by the 
NLSE: 

^) + ^^ + S^)W«)lV«)-0 (66) 
together with the normalisation condition (|43|) . Its solution is given by 



with the eigenvalue 



where 



w 2 cosh(K M x) 



A,, " " 



2m M 



37iiu(l — f 2 ) 



(69) 



Thus, from (|55|) we find that 



h 2 K 2 

= E^kaJ - Vk Cfl - —a-. (70) 
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= \/lf y, t/j r- (71) 

V o eosh/c^n + ±j — Vt — zq) 

This excitation is spatially distributed between the chains with the probability components given 
by: 

Pj,n(t) = l<Pl(™+ J ~-Vt-z ). (72) 



Clearly, Pj = ^ n Pj : n = 1/3. For the totally symmetric soliton, \i = 0, the chains are excited 
with the same phase, while for the other two cases, [i = ±2tt/3 and the excitations in the spines 
have the phase shifts ±27r/3. 

Note that due to the factor (1 — v 2 ) in (|65[) and ()69|) we see that the soliton velocity V cannot 
exceed the sound velocity v a . However, there is also a further restriction on the soliton velocity 
which follows from Q53j) . Unlike for the parabolic law, the energy dispersion in an exciton band 
shows that dE(k)/dk has a maximum value. Therefore, (|53|) has a solution only when V does 
not exceed the maximum exciton group velocity V g = (l/h)(dE(k)/dk) max and so, the top speed 
of the solitons is determined by the lowest of v a and V g . For example, in a simple chain with 
E(k) = —2Jcosk, V g = 2J/h, and with the parameters of the a-helix, v a > V g . 

Below we will consider solitons at low velocities. In this case, from (jEHJ), we have 

k c » = k^ + ^V. (73) 

Here kn determines the bottom of the fi-th exciton band and is an effective exciton mass 
near the band bottom. At low velocities the total energy = 7i of the soliton state is given by 



£»{V) = £^) + -M,V 2 . (74) 



The totally symmetric exciton band has a minimum at &o = and, in the long-wave approx- 
imation, we have 

9ft 2 

£b(0) = E -2J + 2L, m = 2(gJ _ - (75) 
Therefore, the totally symmetric soliton state is characterized by the width parameter 



the energy 



and the mass 



3x 2 (7R , 

K ° = — 7777 Fv ( 76 ^ 

w(9J — L) 



e ( ) = E -2J + 2L- 3vj2( £_ Ly (77) 

8y 4 

Mo = mo + 3vW(9j-Ly (78) 



The other two of the three soliton states are formed by excitons from the other two bands, 
fj, = ±27r/3. Due to the helical symmetry, the bottoms of these bands are determined by the non- 
zero wavenumbers k± = ±kd- For an a-helix the parameter kd is small and can be determined 
in the long-wave approximation as 

9L 

k d = . (79) 

V3(18J + L) V ; 
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E ± (0) = E 1= E -2J-L- wrrTry m±=mTI = mi , (80) 

Therefore, these soliton states are characterized by the width parameter 

6 X 2 , . 

Kl "u;(18J + L)(l- V 2 )' 18 j 

the total energy at rest 

qt2 9 4 

^(0) = E -2J-L- — — -, (82) 

u ; 2(18J + L) 3w 2 (18J + L)' v ; 

and by the soliton mass 

16y 4 

Mi = ™i+ one, , n 2 2 - ( 83 ) 
3(18J + L)w^v^ 

Note that the energies of these three solitons are split from the bottoms of the corresponding 
energy bands. We should add that the solutions of these soliton states were also found in [P5] . 

The energy levels of the last two solitons are degenerate. However, according to the Jan- Teller 
theorem, this degeneracy can be broken by the distortions of the chains and a hybridization of 
these two states can take place. Below we consider such a case when i.e., ip± ^ and (po = 0. 

In this case we find from (|57|) that ip± are determined by the system of equations: 

A^ + (*) + ~ ~ e-^' k ->V-{x)^{x) = 0, (84) 



In this case, the components Q± of the deformation of the o-helix are also non-zero: 



(85) 



Substituting this into (|58|) we find that, for small velocities, 

v + = -f^ e ~ 4(fc+ ~ fc ~M(*)M*), v_ = v;. (87) 

The deformational potential Vo of the totally symetric distortions is given by 

Vo(*) = -f^(lM*)l 2 + IM*)l 2 )- 



Thus, equations (|51)) and give us a system of NLSEs: 
and, equivalently, for if-. 
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V± = -jf 9± ^ ( 9 °) 
where 0± are arbitrary phases and </?2 satisfies the NLSE and is, therefore, given by (|67|) with 

K2 = M0Trj - (91) 

The total energy of this soliton state at V = is 

Q7-2 o 4 

£,(0) =E „-2J-L-^ JTIy - 2m2(1 ^ + L) , (92) 



and the soliton mass is 

22x 4 
w 2 v%(18J + L) 



M 2 =m l + 2 % ; FT . (93) 



Representing the energies of the other two solitons (j82j) in the form £i(0) = — A with 
Ef, = E±(k c i) being the corresponding bottom of the energy band and 

2 X 4 

A = 3^ 2 (18J + L)' (94) 

we can write 

£ 2 (0)=E b (0)-^A = £ 1 (0)-^A. (95) 
Thus, we see that the latter hybrid soliton has the lowest energy. 

The distribution of the excitation amongst the chains is given by the probability amplitude: 

V'j.nC*) = y 3 e 1 cos I k d {n + j/3) + + —j \ (p 2 {n + -, t). (96) 

Therefore, 

(fc) K2 l + cos(2Mn + J 73) + (^-^)-f J ) 

3 cosh (re + | — V t — zo) 

Next, we consider the probability distribution of the excitation summed over all the spines 
of the helix: 

P n (t) = J2 p jA t ) = \Mn,t)\ 2 (^-0=cos(2k d n + e + -e^ ^\^ 2 {n,t)\ 2 (98) 
for small k d , and the total probability of the excitation localisation on a given spine: 



m =Y, p ut)\ 2 = \ 



Trkd ( 2tt 

r- cos 2k r jVt ■ 

K2 S i n h 2& V 3 ' 



(99) 



We see from (|99j) that the probability of the excitation localisation on a given spine is an 
oscillatory function of time with the period of oscillations given by 

T = W- < 100 » 
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period of oscillations that is determined by the soliton velocity and the quasimomentum value 
corresponding to the bottom of the energy band (|79|). These oscillations get mixed up with the 
oscillations that arise from the influence of the lattice discretness on the soliton dynamics which 
leads to the appearance of the Peierls-Nabarro potential. The period of these latter oscillations 
is also determined by the soliton velocity, , as is shown in ^U] . 



5 Numerical modeling 

For the numerical calculations we consider an a-helical system of length N = 150 with periodic 
boundary conditions: 

fj,n+N = fj,n (101) 

or, equivalently, 

fl+3N = fl, (102) 

where / stands for tp or (5 and the index I ennumerates the sites along the polypeptide chain 
and {j, n} denotes the site number n in the j-th hydrogen bound spine (j = 1, 2, 3). 

It is more convenient, for the numerical simulations, to use the physically more relevant 
site representation for the VPj n variables and to use u^ n for the displacements of PGs from the 
positions of their equilibrium. Here the average displacements of PGs in the state ()40|) 

and are related to /?„ )9 by the unitary transformation (|36|) . In these variables the equations 
(|35 |) - (fi5|) become 



in ^f - = ^O^l.n " e/(*l, n -l + *l,n+l) + £(*3,n-l + ^2,n) + 

X(u 1<n+ i - ui in _i)*i jn , (103) 

i% ~ir = E ^ 2 > n ~ J (^2,n-l + ^2,n+l) + l,n + ^3,n) + 

X(u2,n+1 ~ W2,n-l)*2,n, (104) 
iTl ~ir = E °^ 3 > n ~ J (^3,n-l + *3,n+l) + L(^2,n + *l,n+l) + 

x{u-i, n +i ~ n 3jn _i)^3 in , (105) 

d 2 u- 

M —^T = - W ( 2u i,n - Uj, n -1 ~ U jtn+1 ) + X (| ^n+l | 2 " |*i,n-l| 2 ), j = 1, 2,3. (106) 

For n = and n = N — 1 in the expressions above we take the appropriate values of the functions 
determined by our periodicity conditions. 

In our studies we have adopted the following procedure. We have started off with a rea- 
sonable field configuration and then used it as an initial excitation to determine a stationary 
solution of our system of equations (|103I106|) . Having determined this solution numerically, we 
have kept on modifying it by an adiabatical increase of the wave-vector (thus increasing the 
velocity of the soliton) , and have found for each fixed value of the wave- vector the correspond- 
ing stationary solution describing a soliton which propagates along the helix with an increasing 
non-zero velocity, determined by the gradually increasing values of the carrying wave-vector. 
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for a simple chain, because there are three types of solutions, corresponding to the different sym- 
metries. Studying similar problem for a simple linear chain, we had two equivalent approaches 
of deriving a stationary solution at zero velocity (see, e.g., JH]). Namely, we could start with 
the system of stationary equations and find the solution by minimizing the energy using some 
standard procedures. Another approach would use the non-stationary equations which include 
some dissipation of the energy in the lattice sub-system. Starting with an arbitrary localized 
initial configuration of an excitation, we would find some time later a stationary solution at zero 
velocity. Then this configuration would be modified by adding a small carrying wave-vector and 
would be used as a starting initial condition for the next set of calculations of the system of 
equations without any dissipation. This would result in a solution moving with a small non-zero 
velocity. Repeating this procedure further, we would increase the soliton velocity adiabatically 
till the velocity reaches the maximum value corresponding to the chosen parameters of the chain. 
This last approach, of using an arbitrary initial configuration, in the case of a helical structure 
probably cannot describe all possible solutions, since it would always lead to the solution of the 
lowest energy. 



The energy expression can be obtained from (j23U21ll2"l)) . In the site representation the total 
energy is given by 

Etot = E e + E v + Ei n t, (107) 

where 



N-l 3 

E e = E E(^l*l?,„ -J(**„%,n-l + **«-l%n)) + 
n=0 j=l 

+ £(*i,n*3,n-l + ^3,n-l^l,n + *2,n*l,n + ^!,n^2,n + 
+ *3,n*2,n + *2,n*3,n)l , 



E r 



3 N-l 

EE 

j=l n =0 



M ( duj n \ 1 2 



3 N-l 

Eint = E E X(Uj,n+l -Uj,n-l)|*||, n - 
j=l n=0 



(108) 
(109) 

(110) 



In our simulations we have taken the numerical values of the parameters from ^2]: i.e. 
L = 12.4cm" 1 , J = 7.8cm" 1 , X = .34- 10" 10 N, w = 19.5 N/m and JWJw = l/u a = 0.99- 10" 13 
s. These parameter values correspond to the Amid-I excitations in a-helices |151 1201 IT2l 15], 



In our numerical studies we have also followed the conventions of Scott |T2] and so, like him, 
we have used units in which the energy is measured in units of hv a , time in units of v~ l and 
length in units of 10" 11 m. In this case the dimensionless computer values of the parameters are 



J 



at X x 10 n m 

Jcomp — ~z — 0.145, L com p — 0.231, Xcomp — z — 0.318, 

hv a * hv a 



w, 



comp 



ro x 10 22 m 2 



1.825. 



(HI) 



The results of the numerical simulations are described below. 
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= 1. 



(112) 



with all other values of $ 3 - n = and putting Uj n = 0. We have then added an extra absorptive 
term into the equation for Uj^ n and performed the simulations until we have reached a stationary 
state solution. The obtained solution described a well defined solitonic state. Its energy was 
around -0.55067. In Fig. 2 we present the plots of the and u fields. 




2. a 2.b 

Figure 2: Stationary excitation function \& (a) and the derivative of the PGs displacement i.e. 
Uj+i,n — Ujn (b) for the alpha-helix. 

We see that this self-trapped state has an inner structure. While the total (summed over 
all three spines) distribution of the excitation has a single-hump pattern, the distributions in 
individual spines are modulated in the manner of solutions 1)9 7|) . The same feature can also be 
seen in Fig. 6 of ^2j. Thus, we can conclude that our numerical solution, as well as the solution 
discussed in [12], describes the lowest energy of the hybrid solitons. This view is confirmed also 
by the numerical estimate of the soliton energy (|92|) . Thus taking our numerical values (jlllj) . 
we get £2(0) — Eq = —0.55062 in units of hu a which coincides with the value determined in our 
numerical simulations. 

Having found stationary solutions, we then changed the functions as follows 

9 j>n - ^ e inAk - ^ + {u hn - %>n _!) sin(Afc) (113) 

leaving Uj n unchanged. This had the effect of giving a small speed to the soliton, and the 
distortion of the chain. 

We then performed the simulation over a short period of time. During this time the soliton 
has been moving and the small disturbance introduced by the nonperfect transfer of momentum 
to the system has spread itself over the lattice. 

We then repeated the whole process several times thus slowly increasing the total k (in 
practice we put AA;=0.1). After every step we evaluated the resultant speed of the soliton. Of 
course, the whole process suffered by the introduced disturbances; thus gradually it has become 
more and more difficult to determine this speed. However, we have found that each addition of 
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Figure 3: Speed of the hybrid soliton as a function of k. 

momentum increased this speed by a decreasing amount suggesting that there is a maximum 
speed that the soliton can attain. 

In Fig. 3 we present a plot of the resultant speed as a function of the total k (i.e., the sum of 
all AA;). We note that the maximum speed appears to be around 0.21. To check that this limit 
is not an artifact of our procedures, we have performed further simulations in which we modified 
the steps A A; or eliminated the modification of ^jf. We have also performed some simulations 
with absorption: the configurations were alternatively boosted and then evolved in time but 
with a small absorption parameter added to the equations. These extra terms absorbed some 
of the ripples while the boosts were effectively accelerating the solitons. All these procedures 
produced similar results and we have never managed to get the solitons move faster than with 
v ~ 0.21. The absorptions did decrease the deformations of the a- helix but they did also reduce 
the velocity of the soliton; hence we do believe that the solitons cannot have larger velocity and 
that this maximum speed is determined by the maximum allowed group velocity of the excitons. 

In Fig. 4 we present the plot of the solutions of for E±{k) with the values of J and L 
given in (|111|) (we recall that v = V/v a ). From Fig. 4 we see that, indeed, the composite soliton 
cannot have its velocity larger than the maximum group velocity for one of its two components, 
and for our parameters this velocity is about 0.21. At wavenumber k cr , which corresponds to the 
maximum group velocity, d 2 E ^(k) / dk 2 = and at k > k cr the balance between the nonlinearity 
and the dispersion breakes down for one of the components and this leads to the decay of the 
soliton. 

The complex (modulated many-hump) and composite (three-spine distributed) structure of 
the soliton manifests itself distinctly when the soliton is moving and the interspine oscillations 
take place (|99[) . This is seen very clearly in the oscillations of the probability distribution 
amplitude for each spine which is shown in Fig. 5. This phenomenon was already noted by 
Scott in • According to (JlOOj) , the frequency of these oscillations is determined by kd and by 
the soliton velocity. It follows from ()79|) that the bottom of the band is attained at = 0.42. 

We have also looked at the other two solitons and tried to make them move. As has already 
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Figure 4: The excitation velocity (|53f) for J = 7.8cm 1 and L = 12.4cm. 




Figure 5: Oscillation frequency of the hybrid soliton as a function of its speed V. 
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leads to the solution of the lowest energy. But when we take as an initial condition \Pj n in the 
form (|71|) at V = 0, we have obtained, as a result of the calculations, stationary solutions and 
these solutions were very close to those derived in the continuum approximation. The energies 
of these stationary solutions were +0.171 for /i = and —0.54972 for [i = ±2tt/3 which, again, 
coincide with the values S ^ (0) — Eq estimated from (|77)l and (|%2"|) . For these states the probability 
distribution in individual spines |^|| n shows a one- hump pattern without any modulation (see 
fig 6). This differentiates these states from the lowest energy composite soliton. 
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6. a 6.b 
Figure 6: Electron fields for the other solitons - (71) with V = 0. 

We have tried to make these states move. Unfortunately, the perturbations introduced by 
the discreteness of the lattice and by the inexactness of our procedure led to their instability. 
This showed itself in the system evolving into the lowest energy (totally symmetric) soliton. 



6 Conditions for the applicability of the adiabatic approxima- 
tion 



Having found the three types of solutions described in the previous sections a question then arises 
about the conditions of the applicability of the adiabatic approximation in such a three-spine 
model. 

The Hamiltonian (|25|) that describes the states of quasiparticles which interact with phonons, 
does not have an exact solution. The adiabatic approximation describes the soliton-like states of 
large polarons when the autolocalization, within the region of several lattice sites, takes place. 
This is one of the three possible approximations which allow us to represent the Hamiltonian 
(|25|) as a sum of two terms: the main part, Hq, and the term, H\, which can be considered 
as a small correction, and, therefore, for which the perturbation theory can be developed. The 
other two approximations correspond to the almost free quasiparticles and to small polarons. 
The realization of one or another of these three regimes depends on the relation between the 
parameters of the system. In general the problem can be investigated in the framework of the 
variational approach |211 1221 123j . The ground state diagram for a simple chain with one exciton 
band and one phonon mode was presented in This diagram showed the range of values of 
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for which one or the other regime was realized. 



As it has often been mentioned, various properties of the Davydov solitons in a-helical 
proteins have been analysed using a single chain model. Although such a model gives good 
qualitative and sometimes also good quantitative J2j properties of Davydov solitons, the ground 
state diagram |211 I23| shows that the parameters of the a-helix applied to a one-chain model, 
correspond to the state far from the region where the soliton ground states are realized. This 
is one of the reasons why the estimates by H. Bolterauer (Hj and J.W. Schweitzer and J. P. 
Cottingham [5] of the Davydov soliton life-time, obtained within a different approach but still 
based on the one-chain model, give very small values. 



Here we return to this problem and we assess the conditions of the aplicability of the adiabatic 
approximation for the a-helix basing our discussion, for simplicity, on the solutions describing 
solitons at rest. Applying a unitary transformation, the Hamiltonian (|25f) takes the form H = 
H a( i + H na where H a d is diagonal in the new represantation and describes the adiabatic states 
of the exciton (electron)-phonon system. The term H na is an operator of nonadiabaticity which 
describes phonon-induced transitions between adiabatic states. Such a transformation was used 
in Uni and, based on it, a method of partial diagonalization was further developed in [2*Hl 

nam. 

The partial diagonalization shows clearly that the state-vector (|1U|) is an eigenstate of H aa - 
with the eigenenergy £ = hQ + W (here W is the energy of lattice deformation) provided that the 
functions ij)^ and j3 u>q are stationary solutions of equations and (|46|). i.e. H a d\ipo) = £s\i>o}- 
The virtual excited adiabatic states, for a given chain deformation (|52[) can be found from the 
linear equation (|45|). 



If H na is small it is possible to construct the perturbation series 

|^} = |^o) + |Vi} + - 

where \iJjq) = \s) is the wavevector (|40|) of the soliton state in the zero order of the adiabatic 
approximation and |V>i) is the i-th. correction due to H na . According to the general theory of 
perturbations the first correction is given by 

= -®H na \ip ) (114) 

where we have defined 

« n a d — o s 

and 

Q = 1- |^o)(*0o| = \ a )( a \- 

Note that here \a) ennumerates all adiabatic terms of H a d, H a d\a) = £ a \ot). For the conver- 
gence of the perturbation series should be proportional to A 8 with A being a small parameter. 

The square of the norm of vector \tp) is 

<V#) = (ipotyo) + m^i) + ... = 1 + 0(A 2 ). 
Therefore, the applicability of the adiabatic approximation is guaranteed provided that 

(Vi|Vi> = A 2 « 1. (115) 
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= {ipo\H na ^H na \ip ) = ^2(ipo\H na \a)(a\Q\a)(a\H na \ip ) (116) 



a 



= ^2 ^2fa(ipo\Hna\a)(a\H na \'ipo). (117) 
Here we have taken into account the fact that the operator ^ is diagonal and we have defined 

f a = -(a\^-\a) (118) 
with A being the energy gap between the solitonic energy level and the lowest excited one. In 

dm 

cr = ^(a|^^|a) (119) 

. CL 

so that J2 a ^s fa = 1- In (|117|) the summation does not include a = s because the diagonal 
matrix elements of the nonadiabaticity operator vanish. Next we observe that 

E f a (MHna\a){a\H na \i; ) < E<V>o|#n» (a\H na \^) = (^\H 2 na \^). (120) 

a a 

Moreover, it is easy to see that 

{MHLWo) = (V'ol^lV'o) - <V>o|^o) 2 = A£ 2 . (121) 



Thus we can derive some estimates without the partial diagonalization of Hamiltonian (|25j) . 
In particular, we can calculate AE 2 using the soliton wavefunction ()4U|) in the zero order adia- 
batic approximation. This way we can estimate the soliton life-time in one chain and we get the 
same result as that obtained by J.W. Schweitzer and J. P. Cottingham [Hj who calculated AE 2 
performing the partial diagonalization, and by Bolterauer |Hj who calculated AE 2 . 

So, the condition of the applicability can be writen as 

<lhhfc><^^<l. (122) 

Calculation of AE 2 gives us 

AE2 4( i_ ^HW), (123) 

where 

_ 4hx 2 sin 2 q _ 16frv a x 2 

3MNuj q 3vru> ' 1 ' 



Taking into account (|52[) we can rewrite (|123j) in the form 

2/^^ 2 sin 2 q I 

ae2 = E 3MN i - 1 E vv, fc vw +9 i 2 

which corresponds to Bolterauer's [0] expression after the transformation to the site representa- 
tion. 
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where 



rN/2 vn 
J-N/2 

Therefore 



J**\ Vll (x)\'dx=- (126) 

at/2 ^ 2k m smh 



For the hybrid soliton we have 



2ixsing / f N / 2 j qx , ,1,2, , f N ^ 2 ,„.,■ „ , .12, 



Qofa) = —7== / e^|^ + (x)|^x+ / eV\ V -(x)\'dx (128) 

W g V3MA^ \J-N/2 J-N/2 J 

and 

/•AT/2 

Q+(g) = Q*U-q) = / e 4 ^)^;(x)^_(x))dx, (129) 

J-N/2 

2 8hv a X 2 n 7.2 2 7T 2 

AE 2 = — (1 j k 2 — — k.2 sm kd cos kdj- (130) 



In the one-dimensional case, the soliton level (121) is a single bound level in the lattice 
deformation potential. Excited adiabatic states belong to the quasi-continuum spectrum with 

eigenenergy A(fc) = which is separated from the soliton level by a gap A = 2 ^ . Therefore, 
we can estimate a 1)119(1 as 

4 

° = N \ p + fc2)2 = — ' ( 131 ) 

where n = 1 for the totally symmetric soliton, and n = 2 for the hybrid soliton since, in this 
case, there are two degenerate bands. 

Finally, we have the conditions for the realization of the soliton-like states: 

(132) 



. 2C w^hv a {9J - L) 
= \l < 1 

for the total symmetric soliton, and 



7; x 2 



Xi _ f^y WOTTX) < 1 (133) 



for the composite soliton, respectively. Here 

1. 



and 



Co — 1 T^o 

7T^ 



C\ = 1 - ^k d K 2 - ^1^2 



where we have assumed that kd <C vr and % < 1. 



23 



llVJliV/ \l 11U U U11>J \j UlVil VV 1U11 UllVj V> v_y X L V_l 1 LI1U11 VV 111V/U V* 1 t. 111 lk_/ V/ UMUUIlllVa 111 Ui.1V> 

partial diagonalization scheme for the one chain model |271 I26| . This condition indicates that 
solitons can exist in soft enough chains and at a strong enough exciton (electron)-phonon cou- 
pling they are stable against quantum fluctuations. The relation (|132[) is the inverse of the 
condition for the weak coupling regime. 

The numerical values of the parameters for the a-helix are: J = 1.55 • 10 -22 Joule, L = 
2.46 • 1CT 22 Joule, x = 35 - 62 • 10~ 12 N, w H = 13 - 19.5 N/m. We can take l/u a = 1CT 13 s. 
For these parameters we get Ao = 2.3 — 11 for the total symmetric soliton, which corresponds 
to the one-band model. Therefore, in this case, the adiabatic approximation is not valid, and, 
consequently, the soliton is destroyed by quantum fluctuations. The corresponding estimates for 
the composite soliton give the value Ai rj 1. For instance, for x = 62 • 10~ 12 N and wh = 14.6 
N/m we get A2 = 0.87, and, therefore, the perturbation series converges. It is worth adding here 
that the larger values of the coupling and the condition that the chains are softer strengthen 
the condition for this type of soliton solution to exist. 



7 Conclusions 



As we have mentioned above, the main aims of this paper were the study of the soliton states 
in a-helical proteins taking into account their helicity structure and the understanding of the 
origin of the inter-spine oscillations observed numerically by Scott in ^2]- The soliton states 
in a-helical proteins are described by a system of nonlinear equations (jl03lll06j) . In our study 
we have restricted the Hamiltonian of amide excitations to two main terms, namely, those that 
describe the intra- and inter- spine interactions, while Scott considered ten additional terms of 
long-range resonance interactions. Our results broadly reproduce the results of Scott. However, 
there are also some differences, which we summarize below. 

The velocity of the soliton propagation in the numerical calculations carried out by Scott 
in |12| . was reported as V = |f a , while our results give the maximum value V = 0.21u a . This 
is due to the fact that in ^2] further terms of the resonance interaction of Amid-I vibrations 
were included, which increase the width of the exciton bands and, therefore, increase the exciton 
group velocity. The additional terms in the Hamiltonian also change the corresponding value of 
kd, but, probably, this change is less significant than the change of the maximum group velocity. 
Nevertheless, our formula (|10fl)l of the period of oscillations for the values v a = lO 13 ,^ 1 and 
kd = 0.42 for the a-helix at V = |f a , gives T = 1.995 • 10 _12 s, which practically coincides with 
the value obtained by Scott, T comp = 2 • 10~ 12 s. 

Our analytical study and the numerical simulations elucidate the conditions for the existence 
of various types of soliton solutions: single-band and mixed two-band solitons. The entangled 
two-band (hybrid) solitons break spontaneously the translational and rotational symmetries, 
and possess the lowest energy. Single-band solutions break only the translational symmetry and 
preserve the rotational symmetry. Single-band solitons turn out to be dynamically unstable: 
once initially formed, they decay rapidly while propagating. There are two main reasons for 
this, which arise from the helical structure of the system, namely, the absence of the forbidden 
gap in the energy spectrum (see Fig. 1) and the Umklamp processes. The absence of the energy 
gap allows the transition to the lowest energy state via the interactions with low-energy phonons. 
The helical symmetry leads to the relation if) IJi (k±2ir) = Vv±i(^)' i- e -' the mixing of single-band 
states takes place, and, as a result, the single-band solutions decay. This is the reason why 
given any initial condition the excitation localises into the state which corresponds to the lowest 
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managed to observe such solitons due to the very special choices of the initial excitations which 
were very close to the expression for the stationary single-band soliton at rest. 

It is also worth comparing this type of solution in a helical system with those in a three- 
spine model without helical symmetry. In the latter case there is a forbidden gap in the energy 
spectrum between the two degenerate bands and the third band. As a result, the initial excitation 
with the energy above the forbidden gap, is self-trapped in a single-band soliton state. The 
totally symmetric soliton predicted analytically in [HI EH was observed numerically in jlOl lllj . 
Such a soliton in a chain without helical symmetry can be destroyed only if a large amount of 
energy is supplied to the system. Therefore, these single-band localised solutions are much more 
stable dynamically than single-band solitons in chains with helical structures. This constitutes 
a qualitative difference between the three-chain system with an helical symmetry and the one 
without it. 

The important question about the existence of Davydov solitons in a-helical proteins remains 
open. Unlike the case of conducting polymers, for which there is reliable experimental evidence 
for the soliton (large polaron and bipolaron) existence, such data are absent for polypeptides. 
The answer to this question is related to the applicability of the adiabatic approximation, which 
is determined by the numerical values of the parameters of a given system. Solitons can exist 
in protein macromolecules provided their parameters satisfy the condition of the adiabatic ap- 
proximation. Note, e.g., that the spring constant for the hydrogen bond Wjj was determined in 
|20| to be 21 N/m. Scott j!2j . who takes into account that the hydrogen bonds in the a-helix 
are 22° oriented, uses the value 19.5 N/m. But as it has been shown above, the effective value 
is w = ^ 2 wh where 7 is determined in (|19j) . For the parameters of the a-helix a = 5AA and 
r = 1.7A we get 7 = 0.9 and, therefore, w = 17.05 N/m. Thus, the geometrical factor helps to 
satisfy the condition for the exsistence of a soliton. 

As we have seen above, the generally accepted parameters for Amid-I excitation do not 
favour the existence of single-band solitons. On the other hand, they are proper for the existence 
of the entangled soliton states, although the nonadiabatic corrections are also important and 
ought to be taken into account. Thus, the one-chain model can give good qualitative results, 
but conclusions concerning the existence and stability of soliton states, based on numerical 
calculations within such an oversimplified model, may not always be correct. Of course, our 
estimates are relatively rough, and the method of partial diagonalization of the Hamiltonian 
would provide better results. Its generalization to systems involving three-chain macromolecules 
can face the problem of the applicability of the long-wave approximation. In such cases the 
partial diagonalization method developed for discrete models by Clogston et al. [25] may turn out 
to be useful. Moreover, the variational methods can give better results (see, e.g., [29 l 122} l2T | 123] ) 
for the crossover states when the perturbation scheme parameter is not very small. 
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